Vortex phase diagram in rotating two-component Bose-Einstein condensates 



o 
o 

(N 

o 



> 
in 

m 
O 
m 
O 



c 

o 
o 



Kenichi Kasamatsu 1 , Makoto Tsubota 1 , and Masahito Ueda 2 
1 Department of Physics, Osaka City University, Sumiyoshi-Ku, Osaka 558-8585, Japan 
2 Department of Physics, Tokyo Institute of Technology, Meguro-ku, Tokyo 152-8551, Japan 

(Dated: February 2, 2008) 

We investigate the structure of vortex states in rotating two-component Bose-Einstein condensates 
with equal intracomponent but varying intercomponent coupling constants. A phase diagram in 
the intercomponent-coupling versus rot at ion- frequency plane reveals rich equilibrium structures of 
vortex states. As the ratio of intercomponent to intracomponent couplings increases, the interlocked 
vortex lattices undergo phase transitions from triangular to square, to double-core lattices, and 
eventually develop interwoven "serpentine" vortex sheets with each component made up of chains 
of singly quantized vortices. 
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A system described by a multicomponent order param- 
eters can excite various exotic topological defects. There 
have been great interest in the study of such topological 
defects, which has been explored extensivly in condensed 
matter systems such as anisotropic superfluid 3 He 0. 
This study is also closely related with some problems in 
superconductors of UPt 3 and SroRuC^ which are thought 
to have non-s-wave symmetry |2| and high-energy physics 
and cosmology Dilute atomic Bose-Einstein conden- 
sates (BECs) offer another superior testing ground to 
investigate the physics of exotic topological defects be- 
cause they are free from impurities and well controlled 
by optical techniques. A trapped BEC in a rotating po- 
tential generates quantized vortices which are topological 
defects characteristic of superfluidity, and several groups 
have suceeded in observing ordered vortex lattices [4j. 
While many interesting phenomena have been found in 
single-component BECs -4], a rich variety of static and 
dynamic phenomena are expected to occur in a system 
of rotating two-component BECs consisting, for example, 
of two different hyperfine spin states of atoms 0. This 
system is characterized by three coupling constants de- 
noted by Cn and C22 (for intracomponents) and Cyi (for 
intercomponent) 0BEE1OIJ] 

In this Letter, we report 
equilibrium structures of vortex states in two-component 
BECs with Cxi = C22 = C but varying values of C\i and 
the rotation frequency based on the numerical analysis 
of the Gross-Pitaevskii equation. 

Mueller and Ho studied the vortex lattice structure 
of two-component BECs by assuming the lowest Lan- 
dau level approximation and a perfect lattice . They 
proposed a phase diagram on the lattice structure and 
found that as Cn/C is increased, the system undergoes 
a continuous structural change from triangular to square 
lattices. Kita et al. discussed the lattice structure in 
F = 1 spinor BECs 12] . However, both studies assume 
the fast rotation limit where the rotation frequency is 
very close to the radial trapping frequency, so a vast ex- 
perimentally accessible parameter regime remains to be 
explored. When there is no rotation, it is known that 



we obtain miscible condensates as the ground state for 
C12/C < 1, while for C12/C > 1 the two components 
phase separate 0, 0, E3 • In the presence of rotation, our 
numerical analysis reveals a rich variety of vortex states 
which have eluded analytic treatments. In particular, 
we find that for C12/C > 1 double-core vortex lattices 
or interwoven "serpentine" vortex sheets become stabi- 
lized for a considerable range of rotation frequency; these 
states can therefore be observed in current experimental 
situations. 

The equilibrium solutions of two-component BECs in 
a frame rotating with the angular velocity fi = ilz 
are obtained by solving the time-independent coupled 
Gross-Pitaevskii (GP) equations for the condensate wave 
functions (hi + Y,j=i,2 Cij|*j| 2 )*i = (i = 1,2) 

00,11 HQ, wh ere h\ = -(m 12 /m t )V 2 + VS - (Sl/Q)L Z 
with the reduced mass m\i — mi 7712/(^1 + m^) and 
the harmonic potential Vi = (rriiwf / '4mi2W 2 )r 2 . Here, 
length and energy are measured in units of &ho = 
y^h/2mi20j and Tluj = %(u>i + 0)2)/%, respectively. In 
this work, the numerical calculation is done in the two- 
dimensional x — y space. The wave functions are nor- 
malized as J dxdy\^i(x, y)\ 2 — 1 and, then, the intra- 
and intercomponent coupling constants are written as 



C u = 87r(m 12 /m i )A'' i 



2d, 



and Gi, 



inNf d a 12 (i ? j) 



with the corresponding s-wave scattering lengths an, a 2 2, 
012 (assumed to be positive), and the particle numbers 
N 2d per unit length along the z-axis. 

For simplicity, the number of the parameters is reduced 
by assuming an = 022, mi = 1112 and loi = L02 — u), and 
N'f d — N% d — N 2d . For the typical experimental con- 
ditions (ujj_,uj z ) — 27r(8,5)Hz and N ~ 10 6 , the values 
of the intracomponent coupling constants are chosen as 
Cxi = C22 = C = 2000. In the numerical calculation, the 
equilibrium solutions are found by the norm-preserving 
imaginary time propagation of the time-dependent cou- 
pled GP equations, starting from arbitrary trial wave 
functions without vorticity. The propagation continues 
until the fluctuation in the norm of the wave function 
becomes smaller than 10~ 8 . 
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FIG. 1: Q-8 phase diagram of the vortex states in rotating 
two-component BECs. The classification of the structure is 
made by the symbols, A: triangular lattice, o: square lat- 
tice, <g): stripe or double-core vortex lattice, o: vortex sheet. 
Because of the continuous change from triangular to square, 
their boundary is shown by using both A and o. The plots 
at Q, = 1 show the results of Ref. The shaded region 

shows the analytically obtained region of vortex sheets; the 
solid curve represents b — A p and the dashed curve 2b — Rtf 
(see text). 



With a fixed value C = 2000, we investigate the equi- 
librium solutions for various values of two free parameters 
5 = G\ijC and f2 = f2/£D. The obtained vortex structure 
is summarized in the phase diagram of Fig. ^ Because 
the centrifugal potential overcomes the harmonic trap- 
ping potential for D, > 1, the upper limit of the rota- 
tion frequency is set by Q — 1. Previous studies have 
revealed that two non-rotating components with equal 
mass and equal intracomponent scattering length over- 
lap completely for 6 < 1, and phase-separate for 6 > 1 
0, 0, UJ|- This criterion is reflected on the structure 
of vortex states, as explained below. In general, there 
are many metastable states in rotating two-component 
BECs; they are obtained by the numerical simulation 
starting from different wave functions. Here we discuss 
some characteristic features of vortex structures. 

In the region 5 < 1, two types of the regular vortex 
lattices are obtained as the equilibrium state. For 6 = 
where two components are not interacting, the formu- 
lation is equivalent to that of one component, thereby 
triangular vortex lattices are formed [l^- As 6 increases, 
the positions of vortex cores in one component gradually 
shift from those of the other component and the triangu- 
lar lattices are distorted. Eventually, the vortices in each 
component form a square lattice rather than a triangular 
one. The S dependence on the stable region of a square 
lattice was studied in Ref. 0] in the high rotation limit 
f2 ~ 1. As seen in Fig. ^ however, that stable region 
depends not only on S but also sensitively on Q. An 
increase in rotation frequency also causes the transition 
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FIG. 2: Two-dimensional density profile of the condensates 
with C = 2000 and 8 = 0.7 for (a) fi = 0.6 and (b) fi = 0.75. 
(c) Cross section of (b) along the y = line. The total density 
pr = \^i\ 2 + | 2 1 2 and the Thomas-Fermi density profile are 
also shown. 



from triangular to square lattices; in Figs. [21(a) and (b), 
the two-dimensional density profiles of the condensates 
with S = 0.7 are shown for Q = 0.6 and fl = 0.75. We 
find that the two vortex lattices are interlocked in such 
a manner that a peak in the density of one component is 
located at the density hole of the other, as shown in Fig. 
121(c). As a result, the total density p T = |*i| 2 + |*2| 2 
obeys the Thomas-Fermi distribution applied to the over- 
lapping two-component BECs with solid-body rotation 
p T (r) = 2 v /o7^- ar 2 with a = (1 - tt 2 )/C{l + 5). We 
have confirmed that this feature holds for other parame- 
ter regimes. 

Why is the square lattice stabilized in two-component 
BECs? According to the energy functional of this system, 
two components interact via the intercomponent interac- 
tion Ci2|^'i| 2 | v I / 2| 2 ; the velocity field in one component 
is independent of that of the other. Therefore, such a 
feature is determined only by the density distribution of 
the condensates which minimizes the interaction energy. 
We rewrite the interaction energy in terms of the total 
density pr and the "spin" variable S — | 'I'll 2 — l^l 2 as 
Em = j drC[(l + 5)p% + (1 - <5)S 2 ]/4. A larger 5 makes 
a smoother total density px more favorable in order to 
reduce the first term of Ei nt . This results in the shift of 
the positions of vortex cores of each component. Then, 
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FIG. 3: The density profiles of the condensates \&i and $2 
for Q = 0.7 and 5=1. The profiles (a) and (b) are obtained 
by the numerical simulation starting from different trial wave 
functions. The energy deference between these states is AE ~ 
10" 5 . 



the last term of E- lnt can be interpreted effectively as the 
interaction energy between spins in the Ising model |l4| . 
Spin-up components correspond to the density peaks of 
f ! at the vortex cores of Vf^, and vice versa for spin- 
down components. When the coefficient 1 — 6 is posi- 
tive, the interaction between spins is anti-ferromagnetic, 
which makes a square lattice stabilized (because the tri- 
angular lattice could be frustrated). Anti-ferromagnetic 
nature is made more pronounced as the interlocking of 
vortex lattices becomes stronger with 6 and the number 
of vortices increase with Cl. 

As 6 exceeds unity, the system enters a ferromagnetic 
phase. Then, the condensates undergo phase separation 
to spontaneously form domains having the same spin 
component |l4| . Concurrently, vortices of the same com- 
ponent begin to overlap at 8 ~ 1. Figure shows typical 
solutions at 6 = 1. In Fig. Efa), vortices of each compo- 
nent are overlapped in lines, and each condensate den- 
sity forms a stripe pattern . While this structure can 
also be derived if one assumes a perfect lattice ^jji we 
find that its energy is nearly degenerate with that of Fig. 
E^b). For the same parameters, we can obtain another 
equilibrium state called "double-core vortex lattice" in 
Fig. ^ where a vortex lattice of component 2 is made 
by pairs of vortices with the same circulation; vortices 
in component 1 surround those pairs. Therefore, vari- 
ous metastable structures will appear in this parameter 
region. 

For strongly phase-separated region 6 > 1 , the density 
peaks of the same spin component, at which the other- 
component vortices are located, merge further, resulting 
in the formation of vortex sheets; a typical example is 
shown in Fig. ^a). Singly quantized vortices line up in 
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FIG. 4: (a) The density profiles of the condensates for = 
0.75 and 5 = 1.1. The vortex sheets are made up of chains 
of singly quantized vortices whose positions are marked by x . 
(b) The density profile of ^1 and ^2 in the radial component, 
and the corresponding velocity profile. The velocity of solid- 
body rotation is shown by the dotted line, (c) Schematic 
illustration of the model of the vortex sheet state. 



sheets, and the sheets of component 1 and 2 are inter- 
woven alternately from the center to the outward. Fig- 
ure Hfb) shows that each superfluid velocity Vi (i = 1, 2) 
jumps at the vortex sheet, following approximately the 
velocity of solid-body rotation. In the region S > 1, even 
though the value of 6 is changed, the total density is fixed 
by the Thomas- Fermi distribution with 6=1. 

The stationary vortex sheet has been observed in ro- 
tating superfluid 3 He-A 0], where vortices are bounded 
to a topologically stable domain wall soliton across which 
the unit vector representing the direction of the orbital 
angular momentum of the Cooper-pair faces oppositely. 
As discussed in superfluid 3 He-A, the equilibrium dis- 
tance between the sheets of component 1 and component 
2 is determined by the competition between the surface 
tension a of the domain wall and the kinetic energy of 
the superflow 0|. In order to estimate that distance in 
two-component BECs, we consider a simple model shown 
in Fig. IUc). I 11 this model, each velocity is assumed to be 
constant between the vortex sheets of the same compo- 
nent; the value of Vi increases by 2flb across every sheet. 
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The total density pr — | 'I'll 2 + I ^2! 2 is constant, and the 
domain boundary with the penetration depth A(< b) is 
approximated by the linear profile. We calculate the free 
energy F = E — ttL z in the range < r < 2b, the sheet 
distance b being determined so as to minimize F per unit 
area. Firstly the penetration depth A is determined by 
minimizing the surface tension a of a single domain wall 
with respect to A |inj . which is the sum of the quantum 
pressure energy E qp /2nb » (pr/2A) and the interaction 
energy in the overlap region Ei nt /2irb = (8 — l)Ap 2 n /12; 
we thus obtain A p — y/6/C(8 — ljpr- Then, the surface 

tension is written as a — pi^/ 2 '\/C(8 — l)/6. Secondly, 
the flow energy per area (in rotating frame) is given 
by (1/4tt6 2 ) / 2i 'd 2 r^ Pl (v t - tlrf /2 = 29p T ^ 2 /768, 
where b -C A p is assumed. Thus, the free energy per unit 
area is written as 

F 29p T & 2 b 2 2a Cpl 

= — — 1 1 — fll 

4tt 2 768 6 2 1 ' 

Minimizing this energy with respect to b, one obtains 
b = (768cr/29^ 2 p T ) 1/3 = (768/29^ 2 A p ) 1 / 3 . By using the 
Thomas-Fermi density at r = with 6 = 1 as the value 
of pr, the sheet spacing b oc (8- 1) 1 / 6 ^ - ti 2 ) 1 / 12 /tt 2 / 3 is 
consistent with that of numerical solutions; for example, 
for parameters used in Fig. 0] we obtain b = 3.09&ho- 
The vortex sheet is expected in the region b > A p and 
2b < i?TF = \/2/V™i as shown in the shaded region 
in Fig. H When 2b > Rtf, i.e., b becomes comparable 
with the condensate size, the clear structure of sheets 
vanishes. 

Finally, we comment on what structure could be ob- 
served in actual experimental conditions. The two- 
component BECs realized in JILA is a mixture of the 
states 1 1, — 1) and |2, 1) of 87 Rb @. This mixture has 
the scattering lengths which have the same order, an : 
a 2 2 : a 12 = 1 : 0.94 : 0.97, i.e., 8 ~ 1 for equal par- 
ticle numbers. Therefore, lattices with partially over- 
lapping vortices of the same component, as shown in 
Fig. |21 are expected to be observed. The MIT group 
has observed the phase separation of a mixture of 23 Na 
BECs with 1 1,0) and |1,1) state [HI, which has 8 > 1 
(an : 022 : 0,12 = 1 : 1.035 : 1.035). In addition, a 
mixture of 41 K and 87 Rb BECs reported recently [l8| 
lies deeply in a phase-separate region. The vortex sheets 
should therefore be observed at high rotation frequencies. 

In conclusion, we reveal a rich phase diagram of vor- 
tex states in rotating two-component BECs, even in the 
restricted parameter space (8, f2). The structure of a vor- 
tex lattice is shown to depend sensitively on the rotation 
frequency In particular, we find new phenomena such 
as double-core vortex lattices and interwoven serpentine 
vortex sheets. Use of different atomic masses and intra- 
component interactions will realize a coexistence system 
of vortices with different vortex-core sizes; such a situa- 
tion may change the lattice structure shown in this Let- 



ter. We plan to study a more detailed phase diagram and 
effects of the internal Josephson coupling on the vortex 
states. 
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